function L = ThetaLEstimatorMSM(EPointEst,AchievementData,knotvec_c,seed,initialguess,Display,bootstrap,PreliminariesOnly,InitialPopMat)
%%%%NOTES: the first element of bootstrap indicates what method will be used to compute estimates, as follows:
%%%%bootstrap(1) == 0 means compute point estimates via pattern search
%%%%bootstrap(1) == 1 means compute bootstrap estimates via pattern search, with large initial MeshSize=4096 and 'GSS'
%%%%bootstrap(1) == 2 means compute bootstrap estimates via pattern search, with large initial MeshSize=32 size and 'GPS'
%%%%bootstrap(1) == 3 means compute bootstrap estimates via pattern search, with large initial MeshSize=4 size and 'GPS'
%%%%bootstrap(1) == 4 is used for post-estimation of bootstrapped Theta_l types (holding fixed the original dataset, but varying common parameter estimates) 
%%%%bootstrap(1) == 10 means compute point estimates via genetic algorithm and MaxGenerations=10
%%%%bootstrap(1) == 11 means compute point estimates via genetic algorithm and MaxGenerations=200
%%%%bootstrap(1) == 12 means compute bootstraps via genetic algorithm and MaxGenerations=10
%%%%bootstrap(1) == 13 means compute bootstraps via genetic algorithm and MaxGenerations=200


% Display = 0 ;
if nargin<7
    bootstrap = 0 ; 
    PreliminariesOnly = 0 ;
    InitialPopMat     = [] ;
elseif nargin <8
    PreliminariesOnly = 0 ;
    InitialPopMat     = [] ;
elseif nargin <9
    InitialPopMat     = [] ;
end
%%
UdomHet       = EPointEst.UdomHet ; 
F_usmoothHet  = EPointEst.F_usmoothHet ; 

cutoffs       = EPointEst.cutoffs ; 
tau0          = EPointEst.tau0 ; 
tau1          = EPointEst.tau1 ; 
phi           = EPointEst.phi ; 

tol = 10^-6;  %%%%This will be used for monotonicity and convexity constraints
if isempty(Display)
    Display = 0 ;
end

if isempty(seed)
    seed = round(pi*1000000) ;
end

student_id_l = AchievementData.id ; 
Q            = AchievementData.Q ; 
T            = AchievementData.T - AchievementData.Tunfinished ; 
ThetaE       = AchievementData.ThetaE ; 
contract     = AchievementData.contract ; 
base         = AchievementData.b ; 
wage         = AchievementData.w ; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  MATRICES OF SIMULATED PRODUCTION SHOCKS:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if bootstrap(1)==0  %%%%This is the case when doing point estimates
    S         = 160 ;
elseif bootstrap(1)==4  %%%%This is the case when doing Type bootstrap estimation (i.e., holding fixed the original dataset, while varying the common structural parameters within their respective sampling distributions).
    S         = 160 ;
else  %%%%Barring the two cases above, this means that we are doing the initial phase of the bootstrap where we re-estimate common structural parameters repeatedly.
    S         = 40 ;
end
Qbar          = 160 ;  %%%%This is the initial number of simulated production units
NActive       = sum(~isnan(ThetaE) & Q>=2) ;  %%%%This is the number of students who completed at least 2 learning tasks
Activeids     = student_id_l(~isnan(ThetaE) & Q>=2) ;
NMarginal1    = sum(~isnan(ThetaE) & Q==1 & contract==1) ;  %%%%This is the number of students in contract group 1 who completed 1 learning task but not two
Marginalids1  = student_id_l(~isnan(ThetaE) & Q==1 & contract==1) ;
NMarginal2    = sum(~isnan(ThetaE) & Q==1 & contract==2) ;  %%%%This is the number of students in contract group 2 who completed 1 learning task but not two
Marginalids2  = student_id_l(~isnan(ThetaE) & Q==1 & contract==2) ;
%%
stopflag      = 0 ;
notenoughflag = zeros(NActive+NMarginal1+NMarginal2,1) ;
%%%%THIS LOOP INCLUDES THE MARGINAL STUDENTS FROM CONTRACTS 1 AND 2...
while stopflag==0 %%%%This loop will increase the number of simulated production units as needed so that we can be able to match all subjects' ThetaL's
    rng(seed) ;
    SimTimeCum = rand(Qbar,S,NActive+NMarginal1+NMarginal2) ;  
    %%%%This will contain the cumulative sum (within columns)rand(Qbar,S,NActive+NMarginal1+NMarginal2 of SimTime which, when scaled by 
    %%%%ThetaE the qth row represents a simulated cumulative work time through q units of success.  
    %%%%NOTE: FOR COMPUTATIONAL SPEED WE SCALE SIMULATED TIMES BY THETA_E HERE, TO AVOID LOTS OF 
    %%%%REDUNDANT FLOPS WITHIN THE OBJECTIVE FUNCTION EVALUATIONS.
    for ii=1:NActive+NMarginal1+NMarginal2
        if ii<=NActive  %%%%This means we are looking at an active student
            tempid = Activeids(ii) ;
        elseif ii>NActive && ii<=NActive+NMarginal1  %%%%This means we are looking at a marginal student in contract group 1
            tempid = Marginalids1(ii-NActive) ;
        else %%%%Otherwise, we are looking at a marginal student in contract group 2
            tempid = Marginalids2(ii-NActive-NMarginal1) ;
        end
        tempE              = ThetaE(student_id_l==tempid) ;
        tempseg            = find(cutoffs<tempE,1,'last') ;
        if isempty(tempseg)
            tempseg = 1 ;
        elseif tempseg==length(cutoffs)
            tempseg = length(cutoffs) - 1 ; 
        end 
        F_utemp        = F_usmoothHet(:,tempseg) ; 
        Udomtemp       = UdomHet(:,tempseg) ; 
        flag           = find(abs(F_utemp(2:end)-F_utemp(1:end-1))<10^-6 | abs(Udomtemp(2:end)-Udomtemp(1:end-1))<10^-6) ; 
        F_utemp(flag)  = [] ; 
        Udomtemp(flag) = [] ; 
        SimTimeCum(:,:,ii) = interp1(F_utemp,Udomtemp,SimTimeCum(:,:,ii),'linear') ; 
        SimTimeCum(:,:,ii) = SimTimeCum(:,:,ii)*tau0 ; 
        if EPointEst.BurninModel==1 
            SimTimeCum(1,:,ii) = SimTimeCum(1,:,ii)*tau1 ; 
        end 
        learning           = ( ((1:Qbar)').^-phi )*ones(1,S) ; 
        SimTimeCum(:,:,ii) = SimTimeCum(:,:,ii).*learning ; 
        SimTimeCum(:,:,ii) = cumsum(SimTimeCum(:,:,ii))*tempE ; 
        notenoughflag(ii)  = 1*( mean(SimTimeCum(Qbar,:,ii)) <= T(student_id_l==tempid) ) ; 
    end 
    if sum(notenoughflag)>0 
        Qbar = Qbar + 10 ; 
    else 
        stopflag = 1 ; 
    end 
end 
%%




%%%%Consistent with our model of the continuation value based on projected future payoffs of the 
%%%%form E[Q|t]x\Pi(Q):
unit1workshare = (tau1*(1+1-2^-phi))/(1+(tau1*(1+1-2^-phi))) ;  %%%%This is the expected fraction of time for the first two units that is accounted for by the first unit. 

pay1 = unique(base(contract==1)) + unique(wage(contract==1))*(2:Qbar)' ; %%%%This is the pay structure under contract 1 
pay2 = unique(base(contract==2)) + unique(wage(contract==2))*(2:Qbar)' ; %%%%This is the pay structure under contract 2 
pay3 = unique(base(contract==3)) + unique(wage(contract==3))*(2:Qbar)' ; %%%%This is the pay structure under contract 3 

pay  = [[0;unit1workshare*pay1(1);pay1] [0;unit1workshare*pay2(1);pay2] [0;unit1workshare*pay3(1);pay3] (0:Qbar)'] ; 
%%%%This matrix contains pay structure under the various contracts.  Columns 1-3 correspond to contracts 1-3, respectively, 
%%%%while column 4 contains the number of units produced for the pay in each row.  For the first unit of output, we model pay 
%%%%as a partial payment of \Pi_i(2): in expectation, unit1wkfrac is the fraction of work that will go into the first payout
%%%%for 2 completed tasks, so we set \Pi_i(1) = unit1wkfrac x \Pi_i(2).
%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% NOW WE COMPUTE B-SPLINE CDFs FOR Q, INCLUDING AN EXTRAPOLATED UPPER TAIL.  THE EXTRAPOLATION IS
%%%% TO DEAL WITH THE SMALL MASS POINT AT FULL OUTPUT (80): FOR THESE PEOPLE WE DON'T KNOW THEIR
%%%% TRUE OUTPUT; JUST THAT IT WOULD HAVE BEEN AT LEAST 80 IF THEY COULD HAVE CONTINUED.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
UpperTailRank = mean(AchievementData.Q(AchievementData.contract==3)>=2);  %%%%This is the quantile rank of the marginal active student in Contract 3.
C3ESampSize1  = round(UpperTailRank*sum(AchievementData.contract==1));  %%%%This is the sample size of the C3Equivalent in Contract group 1.
C3ESampSize2  = round(UpperTailRank*sum(AchievementData.contract==2));  %%%%This is the sample size of the C3Equivalent in Contract group 2.
C3EAddons1    = C3ESampSize1 - round( sum(AchievementData.Q(AchievementData.contract==1)>=2) );   %%%%This is how many marginal achievement levels to add to sample 1 to get C3 equivalence.
C3EAddons2    = C3ESampSize2 - round( sum(AchievementData.Q(AchievementData.contract==2)>=2) );   %%%%This is how many marginal achievement levels to add to sample 2 to get C3 equivalence.


%%%%THESE ARE THE CONDITIONAL WORK CDFs (CONTRACT 3 EQUIVALENT): 
[Fq1,Xq1]     = ecdf([ones(C3EAddons1,1); AchievementData.Q(AchievementData.contract==1 & AchievementData.Q>=2)]);
[Fq2,Xq2]     = ecdf([ones(C3EAddons2,1); AchievementData.Q(AchievementData.contract==2 & AchievementData.Q>=2)]);
[Fq3,Xq3]     = ecdf(AchievementData.Q(AchievementData.contract==3 & AchievementData.Q>=2)); 



K_q           = 10 ;
uppertrunc    = 160 ;  %%%%This is the assumed finite upper truncation point for the quantile
%%%%functions.  I.e., with uppertrunc=160 we assume that no one would have been willing to complete
%%%%more than 200% of the offered workload of 80 in the 10-day period.
if bootstrap(1)==0 || bootstrap(1)==10 || bootstrap(1)==11
    %%%%HERE WE WILL COMPUTE THE INTEGRATED CDF RATIO TO HELP DISCIPLINE THE EXTRAPOLATED
    %%%%UPPER TAIL:
    commondom12   = unique([Xq1;Xq2]) ;
    Fdiff12 = interp1(Xq1(2:end),Fq1(2:end),commondom12,'next') - interp1(Xq2(2:end),Fq2(2:end),commondom12,'next') ;
    Fdiff12 = sum( (commondom12(2:end)-commondom12(1:end-1)).*Fdiff12(1:end-1) ) ;
    commondom23   = unique([Xq2;Xq3]) ;
    Fdiff23 = interp1(Xq2(2:end),Fq2(2:end),commondom23,'next') - interp1(Xq3(2:end),Fq3(2:end),commondom23,'next','extrap') ;
    Fdiff23 = sum( (commondom23(2:end)-commondom23(1:end-1)).*Fdiff23(1:end-1) ) ;
    FdiffRatio = Fdiff23/Fdiff12 ;
    
    
    [Fq,Xq]    = ecdf([ones(C3EAddons1,1); AchievementData.Q(AchievementData.contract==1 & AchievementData.Q>=2);...
        ones(C3EAddons2,1); AchievementData.Q(AchievementData.contract==2 & AchievementData.Q>=2);...
        AchievementData.Q(AchievementData.contract==3 & AchievementData.Q>=2)]) ;
    
    %%%%HERE WE COMPUTE SMOOTHED B-SPLINE CDFs, WITH EXTRAPOLATED UPPER TAILS
    knotvec_q      = Fq(end-1)*linspace(0,1,K_q)' ;
    knotvec_q      = [0; interp1(Fq(2:end),Xq(2:end),knotvec_q(2:end),'linear')] ;
    knotvec_q      = sort([knotvec_q; uppertrunc]) ;
    
    %%%%CONTRACT GROUP 1 CDF:
    C1            = BsplineBasis3(knotvec_q,Xq1(2:end-1)) ; %%%%BASIS FUNCTION VALUES; THE "REGRESSORS" ARE EMPIRICALLY OBSERVED VALUES OF WORK VOLUME
    d1            = Fq1(2:end-1) ; %%%%"OUTCOME VARIABLES"
    %%%%the next two lines enforce monotonicity:
    A1            = [eye(K_q+3-1),zeros(K_q+3-1,1)]-[zeros(K_q+3-1,1),eye(K_q+3-1)] ;
    b1            = (-(10^-6))*ones(K_q+3-1,1) ;
    
    %%%%CONTRACT GROUP 2 CDF:
    C2            = BsplineBasis3(knotvec_q,Xq2(2:end-1)) ; %%%%BASIS FUNCTION VALUES; THE "REGRESSORS" ARE EMPIRICALLY OBSERVED VALUES OF WORK VOLUME
    d2            = Fq2(2:end-1) ; %%%%"OUTCOME VARIABLES"
    %%%%the next two lines enforce monotonicity:
    A2            = [eye(K_q+3-1),zeros(K_q+3-1,1)]-[zeros(K_q+3-1,1),eye(K_q+3-1)] ;
    b2            = (-(10^-6))*ones(K_q+3-1,1) ;
    
    %%%%CONTRACT GROUP 3 CDF:
    C3            = BsplineBasis3(knotvec_q,Xq3(2:end-1)) ; %%%%BASIS FUNCTION VALUES; THE "REGRESSORS" ARE EMPIRICALLY OBSERVED VALUES OF WORK VOLUME
    d3            = Fq3(2:end-1) ; %%%%"OUTCOME VARIABLES"
    %%%%the next two lines enforce monotonicity:
    A3            = [eye(K_q+3-1),zeros(K_q+3-1,1)]-[zeros(K_q+3-1,1),eye(K_q+3-1)] ;
    b3            = (-(10^-6))*ones(K_q+3-1,1) ;
    
    C = [C1, zeros(length(C1(:,1)),2*(K_q+3));...
        zeros(length(C2(:,1)),(K_q+3)), C2, zeros(length(C2(:,1)),(K_q+3));...
        zeros(length(C3(:,1)),2*(K_q+3)), C3] ;
    d = [d1; d2; d3] ;
    A = [A1, zeros(length(A1(:,1)),2*(K_q+3));...
        zeros(length(A2(:,1)),(K_q+3)), A2, zeros(length(A2(:,1)),(K_q+3));...
        zeros(length(A3(:,1)),2*(K_q+3)), A3] ;
    b = [b1; b2; b3] ;
    %%%%CHECKPOINTS below will be used to enforce ordering among the upper CDF tails for the three
    %%%%distributions.
    Ncheckpoints = 500 ;
    checkpoints  = linspace(knotvec_q(1),knotvec_q(end),Ncheckpoints)' ;
    A = [A;...
        -1 zeros(1,3*(K_q+3)-1) ;...                  %%%%This enforces non-negative CDF values
        zeros(1,K_q+3) -1 zeros(1,2*(K_q+3)-1) ;...   %%%%This enforces non-negative CDF values
        zeros(1,2*(K_q+3)) -1 zeros(1,(K_q+3)-1) ;... %%%%This enforces non-negative CDF values
        1 zeros(1,3*(K_q+3)-1) ;...                  %%%%This enforces lower-end CDF values below 10^-2
        zeros(1,K_q+3) 1 zeros(1,2*(K_q+3)-1) ;...   %%%%This enforces lower-end CDF values below 10^-2
        zeros(1,2*(K_q+3)) 1 zeros(1,(K_q+3)-1) ;... %%%%This enforces lower-end CDF values below 10^-2
        zeros(1,K_q+3-1) 1 zeros(1,2*(K_q+3)) ;...    %%%%This enforces CDF values below 1
        zeros(1,2*(K_q+3)-1) 1 zeros(1,(K_q+3)) ;...  %%%%This enforces CDF values below 1
        zeros(1,3*(K_q+3)-1) 1 ;...                   %%%%This enforces CDF values below 1
        zeros(1,K_q+3-1) -1 zeros(1,2*(K_q+3)) ;...    %%%%This enforces upper-end CDF values above 1-10^-2
        zeros(1,2*(K_q+3)-1) -1 zeros(1,(K_q+3)) ;...  %%%%This enforces upper-end CDF values above 1-10^-2
        zeros(1,3*(K_q+3)-1) -1 ;...                   %%%%This enforces upper-end CDF values above 1-10^-2
        -BsplineBasis3(knotvec_q,checkpoints), BsplineBasis3(knotvec_q,checkpoints), zeros(Ncheckpoints,K_q+3);...  %%%%This line forces the CDF under contract 1 to be above the CDF for contract 2
        zeros(Ncheckpoints,K_q+3), -BsplineBasis3(knotvec_q,checkpoints), BsplineBasis3(knotvec_q,checkpoints)];  %%%%This line forces the CDF under contract 2 to be above the CDF for contract 3
    b = [b; 0; 0; 0; 10^(-6); 10^(-6); 10^(-6); 1; 1; 1; -1*(1-10^(-6)); -1*(1-10^(-6)); -1*(1-10^(-6)); -(10^-6)*ones(2*Ncheckpoints,1)] ;
    %%%%In this last section we create a final constraint to enforce that the tails of the quantile
    %%%%functions maintain constant relative magnitudes as is observed in the upper tail of the
    %%%%empirical distribution just below Qcutoff.  Specifically, between q=0.75 and q=0.95 the ratio of
    %%%%contract 2 quantiles to contract 1 quantiles is roughly constant at about 1.57 and the ratio of
    %%%%contract 3 quantiles to contract 2 quantiles is roughly constant at about 1.35.  Therefore, we
    %%%%add two extra terms to the objective function to enforce this relationship above q=0.95.
    Ncheckpoints = 50 ;

    checkpoints  = linspace(80,uppertrunc,Ncheckpoints)' ;
    Basis        = BsplineBasis3(knotvec_q,checkpoints) ;
    rho = 0.1 ;
    fun     = @(beta) sum( (C*beta-d).^2 ) ...
        + rho*sum((  ((Basis*beta((K_q+3)+1:2*(K_q+3)))-(Basis*beta(2*(K_q+3)+1:3*(K_q+3)))) ...  %%%%This enforces constant tail ratio of 1.57 for contracts 1 and 2 as described above
        ./ ((Basis*beta(1:(K_q+3)))-(Basis*beta((K_q+3)+1:2*(K_q+3)))) - FdiffRatio*ones(size(checkpoints)) ).^2)  ;  %%%%This enforces constant tail ratio of 1.57 for contracts 2 and 3 as described above
    betaq0 = initialguessFq(K_q,'CDF') ;
    
    
    options = optimoptions('fmincon','Display','iter-detailed','MaxIterations',10000,'MaxFunctionEvaluations',200000,'ConstraintTolerance',10^-6,...
                           'OptimalityTolerance',10^-6,'StepTolerance',10^-6) ;
    [betaq,fval,exitflag,output]  = fmincon(fun,betaq0,A,b,[],[],[],[],[],options)   %#ok
    %%%%THE SOLVER TYPICALLY DOESN'T GET THE TERMINAL VALUES EXACTLY RIGHT, SO HERE WE TWEAK THE
    %%%%ENDPOINTS TO BE ZERO AND 1, RESPECTIVELY:
    betaq(1) = 0 ; betaq(K_q+3) = 1 ;
    betaq(K_q+3+1) = 0 ; betaq(2*(K_q+3)) = 1 ;
    betaq(2*(K_q+3)+1) = 0 ; betaq(3*(K_q+3)) = 1 ;

  
    betaq1 = betaq(1:K_q+3) ;
    betaq2 = betaq(K_q+3+1:2*(K_q+3)) ;
    betaq3 = betaq(2*(K_q+3)+1:3*(K_q+3)) ;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%  NOW WITH THE EXTRAPOLATED UPPER TAILS WE COMPUTE FREQUENCY TABLES FOR DIFFERENT LEVELS OF
    %%%%  OUTPUT.  THIS IS THE FINAL OBJECT THAT WILL BE PASSED TO THE OBJECTIVE FUNCTION TO DEAL WITH
    %%%%  THE SMALL MASS POINT AT FULL OUTPUT.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    QuantileFunq1 = linspace(knotvec_q(1),knotvec_q(end),1000)';  %%%%This is the quantile function for group 1
    domq1         = BsplineEval3(knotvec_q,betaq1,QuantileFunq1);
    QuantileFunq2 = linspace(knotvec_q(1),knotvec_q(end),1000)';  %%%%This is the quantile function for group 2
    domq2         = BsplineEval3(knotvec_q,betaq2,QuantileFunq2);
    QuantileFunq3 = linspace(knotvec_q(1),knotvec_q(end),1000)';  %%%%This is the quantile function for group 3
    domq3         = BsplineEval3(knotvec_q,betaq3,QuantileFunq3);
    if Display==1
        figure(33)
        hold on
        Fqplot1 = stairs(Xq1(2:end-1),Fq1(2:end-1),'--b','LineWidth',4) ;
        Qqplot1 = plot(QuantileFunq1,domq1,'--b','LineWidth',2) ;
        Fqplot2 = stairs(Xq2(2:end-1),Fq2(2:end-1),'-.r','LineWidth',4) ;
        Qqplot2 = plot(QuantileFunq2,domq2,'-.r','LineWidth',2) ;
        Fqplot3 = stairs(Xq3(2:end-1),Fq3(2:end-1),':k','LineWidth',4) ;
        Qqplot3 = plot(QuantileFunq3,domq3,':k','LineWidth',2) ;
        for ii=1:length(knotvec_q)
            if ii==1
                Knot_uplot = plot([knotvec_q(ii),knotvec_q(ii)],[0,1],'k') ;
            else
                plot([knotvec_q(ii),knotvec_q(ii)],[0,1],'k') ;
            end
        end
    end
    
    %%%% THE FOLLOWING ARE THE CONDITIONAL CDFS, GIVEN WORK ABOVE 80 UNITS:
    tail1  = QuantileFunq1(QuantileFunq1>=Xq1(end-1)) ;
    Ftail1 = (domq1(QuantileFunq1>=Xq1(end-1))-interp1(QuantileFunq1,domq1,Xq1(end-1)))/(1-interp1(QuantileFunq1,domq1,Xq1(end-1))) ;
    tail2  = QuantileFunq2(QuantileFunq2>=Xq2(end-1)) ;
    Ftail2 = (domq2(QuantileFunq2>=Xq2(end-1))-interp1(QuantileFunq2,domq2,Xq2(end-1)))/(1-interp1(QuantileFunq2,domq2,Xq2(end-1))) ;
    tail3  = QuantileFunq3(QuantileFunq3>=Xq3(end-1)) ;
    Ftail3 = (domq3(QuantileFunq3>=Xq3(end-1))-interp1(QuantileFunq3,domq3,Xq3(end-1)))/(1-interp1(QuantileFunq3,domq3,Xq3(end-1))) ;

   
    
elseif bootstrap(1)>=1 
    %%%%We make the boostraps more stable by simply doing linear interpolation instead of the more
    %%%%involved routine for the point estimates.  Specifically, we draw a straight line between
    %%%%the highest value attained by the bootstrap sample (the highest value below 80, that is),
    %%%%and the maximum of the extrapolated quantile function from the point estimates, whose
    %%%%values are stored in the variable bootstrap (with contracts 1, 2, and 3 being in elements
    %%%%2, 3, and 4, respectively).  This allows the upper tail to still vary with the bootstrap
    %%%%sample, while being more stable and requiring far less babysitting to keep the
    %%%%solver from breaking down with odd-shaped samples.
    domq1         = linspace(Fq1(end-1),1,1000)' ; 
    QuantileFunq1 = interp1([Fq1(end-1),1],[Xq1(end-1),bootstrap(2)],domq1,'linear') ; 
    domq2         = linspace(Fq2(end-1),1,1000)' ; 
    QuantileFunq2 = interp1([Fq2(end-1),1],[Xq2(end-1),bootstrap(3)],domq2,'linear') ; 
    domq3         = linspace(Fq3(end-1),1,1000)' ; 
    QuantileFunq3 = interp1([Fq3(end-1),1],[Xq3(end-1),bootstrap(4)],domq3,'linear') ; 
    
    tail1  = QuantileFunq1 ; 
    Ftail1 = domq1 ; 
    tail2  = QuantileFunq2 ; 
    Ftail2 = domq2 ; 
    tail3  = QuantileFunq3 ; 
    Ftail3 = domq3 ; 
end

%%%%TABULATE THE NUMBER OF PEOPLE WHO COMPLETED 80 LEARNING TASKS IN EACH CONTRACT GROUP:
N80_1 = sum(AchievementData.Q==80 & AchievementData.contract==1) ; 
N80_2 = sum(AchievementData.Q==80 & AchievementData.contract==2) ; 
N80_3 = sum(AchievementData.Q==80 & AchievementData.contract==3) ; 


Q80freq_1 = (80:QuantileFunq1(end))' ; 
Q80freq_1 = [Q80freq_1 nan(size(Q80freq_1))] ; 
Q80freq_1(:,2) = round(interp1(tail1,Ftail1,Q80freq_1(:,1),'linear')*N80_1)  ; 
Q80freq_1(find(Q80freq_1(2:end,2)==Q80freq_1(1:end-1,2))+1,:) = [] ; 
Q80freq_1(Q80freq_1(:,2)==0,:) = [] ; 
for ii=1:length(Q80freq_1(:,2))-1
    Q80freq_1(length(Q80freq_1(:,2))-ii+1,2) = Q80freq_1(length(Q80freq_1(:,2))-ii+1,2) - Q80freq_1(length(Q80freq_1(:,2))-ii,2) ; 
end
%%%%As a final step to reduce computational load later on in the objective function evaluations, we
%%%%coarsen the extrapolated upper tail of the empirical CDF so that steps occur in groups of no
%%%%less than 5 observations, but where there are no more steps than Nsteps total.
temp   = [] ;
temp1  = 0 ;
Nsteps = 10 ;
stepsize =  max(5,ceil(sum(Q80freq_1(:,2))/Nsteps)) ;
for ii=1:length(Q80freq_1(:,1))
    temp1 = temp1 + Q80freq_1(ii,2) ;
    if temp1>=stepsize & ii<length(Q80freq_1(:,1)) %#ok
        temp  = [temp; Q80freq_1(ii,1), temp1] ; %#ok
        temp1 = 0 ;
    elseif ii==length(Q80freq_1(:,1))
        temp  = [temp; Q80freq_1(ii,1), temp1] ; %#ok
    end
end
Q80freq_1 = temp ;



Q80freq_2 = (80:QuantileFunq2(end))' ;
Q80freq_2 = [Q80freq_2 nan(size(Q80freq_2))] ;
Q80freq_2(:,2) = round(interp1(tail2,Ftail2,Q80freq_2(:,1),'linear')*N80_2)  ;
Q80freq_2(find(Q80freq_2(2:end,2)==Q80freq_2(1:end-1,2))+1,:) = [] ;
Q80freq_2(Q80freq_2(:,2)==0,:) = [] ;
for ii=1:length(Q80freq_2(:,2))-1
    Q80freq_2(length(Q80freq_2(:,2))-ii+1,2) = Q80freq_2(length(Q80freq_2(:,2))-ii+1,2) - Q80freq_2(length(Q80freq_2(:,2))-ii,2) ;
end
%%%%As a final step to reduce computational load later on in the objective function evaluations, we
%%%%reduce the extrapolated upper tail to clusters of 5:
temp     = [] ;
temp1    = 0 ;
stepsize =  max(5,ceil(sum(Q80freq_2(:,2))/Nsteps)) ;
for ii=1:length(Q80freq_2(:,1))
    temp1 = temp1 + Q80freq_2(ii,2) ;
    if temp1>=stepsize & ii<length(Q80freq_2(:,1)) %#ok
        temp  = [temp; Q80freq_2(ii,1), temp1] ; %#ok
        temp1 = 0 ;
    elseif ii==length(Q80freq_2(:,1))
        temp  = [temp; Q80freq_2(ii,1), temp1] ; %#ok
    end
end
Q80freq_2 = temp ;


Q80freq_3 = (80:QuantileFunq3(end))' ;
Q80freq_3 = [Q80freq_3 nan(size(Q80freq_3))] ;
Q80freq_3(:,2) = round(interp1(tail3,Ftail3,Q80freq_3(:,1),'linear')*N80_3)  ;
Q80freq_3(find(Q80freq_3(2:end,2)==Q80freq_3(1:end-1,2))+1,:) = [] ;
Q80freq_3(Q80freq_3(:,2)==0,:) = [] ;
for ii=1:length(Q80freq_3(:,2))-1
    Q80freq_3(length(Q80freq_3(:,2))-ii+1,2) = Q80freq_3(length(Q80freq_3(:,2))-ii+1,2) - Q80freq_3(length(Q80freq_3(:,2))-ii,2) ;
end
%%%%As a final step to reduce computational load later on in the objective function evaluations, we
%%%%reduce the extrapolated upper tail to clusters of 5:
temp     = [] ;
temp1    = 0 ;
stepsize = max(5,ceil(sum(Q80freq_3(:,2))/Nsteps)) ;
for ii=1:length(Q80freq_3(:,1))
    temp1 = temp1 + Q80freq_3(ii,2) ;
    if temp1>=stepsize & ii<length(Q80freq_3(:,1)) %#ok
        temp  = [temp; Q80freq_3(ii,1), temp1] ; %#ok
        temp1 = 0 ;
    elseif ii==length(Q80freq_3(:,1))
        temp  = [temp; Q80freq_3(ii,1), temp1] ; %#ok
    end
end
Q80freq_3 = temp ;




if Display==1  


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%              FIGURE OS.3:                 %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    figure(33)
    hold on
    sampextrap1 = [] ;
    for ii=1:length(Q80freq_1(:,1))
        sampextrap1 = [sampextrap1; Q80freq_1(ii,1)*ones(Q80freq_1(ii,2),1)] ; %#ok
    end
    sampextrap1 = [sampextrap1; ones(C3EAddons1,1); AchievementData.Q(AchievementData.contract==1 & AchievementData.Q<80 & AchievementData.Q>=2)] ;
    temp        = Xq1(end-1) ;
    [Fq1,Xq1]   = ecdf(sampextrap1);
    stairs(Xq1(Xq1>=temp),Fq1(Xq1>=temp),'--b','LineWidth',1.5)
    sampextrap2 = [] ;
    for ii=1:length(Q80freq_2(:,1))
        sampextrap2 = [sampextrap2; Q80freq_2(ii,1)*ones(Q80freq_2(ii,2),1)] ; %#ok
    end
    sampextrap2 = [sampextrap2; ones(C3EAddons2,1); AchievementData.Q(AchievementData.contract==2 & AchievementData.Q<80 & AchievementData.Q>=2)] ;
    temp        = Xq2(end-1) ;
    [Fq2,Xq2]   = ecdf(sampextrap2);
    stairs(Xq2(Xq2>=temp),Fq2(Xq2>=temp),'-.r','LineWidth',1.5)
    sampextrap3 = [] ;
    for ii=1:length(Q80freq_3(:,1))
        sampextrap3 = [sampextrap3; Q80freq_3(ii,1)*ones(Q80freq_3(ii,2),1)] ; %#ok
    end
    sampextrap3 = [sampextrap3; AchievementData.Q(AchievementData.contract==3 & AchievementData.Q<80 & AchievementData.Q>=2)] ;
    temp        = Xq3(end-1) ;
    [Fq3,Xq3]   = ecdf(sampextrap3);
    stairs(Xq3(Xq3>=temp),Fq3(Xq3>=temp),':k','LineWidth',1.5)

    legend([Fqplot1 Qqplot1 Fqplot2 Qqplot2 Fqplot3 Qqplot3 Knot_uplot],{'G_{a}(a|Contract 1)','Extrapolated Upper Tail, Contract 1',...
        'G_{a}(a|Contract 2)','Extrapolated Upper Tail, Contract 2',...
        'G_{a}(a|Contract 3)','Extrapolated Upper Tail, Contract 3', 'B-Spline Knots'},...
        'Location','SouthEast','FontWeight','bold','FontSize',16)
    xlabel('COMPLETED LEARNING TASKS, A_i','FontSize',24,'FontWeight','bold')
    ylabel('CDFs','FontSize',24,'FontWeight','bold')
    box on
    grid on
end

Q80freq                               = nan(max([length(Q80freq_1(:,1));length(Q80freq_2(:,1));length(Q80freq_3(:,1))]),2,3) ; 
Q80freq(1:length(Q80freq_1(:,1)),:,1) = Q80freq_1 ; 
Q80freq(1:length(Q80freq_2(:,1)),:,2) = Q80freq_2 ; 
Q80freq(1:length(Q80freq_3(:,1)),:,3) = Q80freq_3 ; 
L.domq1         = domq1 ; 
L.QuantileFunq1 = QuantileFunq1 ; 
L.domq2         = domq2 ; 
L.QuantileFunq2 = QuantileFunq2 ; 
L.domq3         = domq3 ; 
L.QuantileFunq3 = QuantileFunq3 ; 
L.Q80freq       = Q80freq ; 

clear sampextrap1 sampextrap2 sampextrap3 temp ii Q80freq_1 Q80freq_2 Q80freq_3 options
clear temp temp1 Nsteps stepsize Xq1 Fq1 Xq2 Fq2 Xq3 Fq3 N80_1 N80_2 N80_3 tail1 tail2 tail3 Ftail1 Ftail2 Ftail3
clear domq1 domq2 domq3 QuantileFunq1 QuantileFunq2 QuantileFunq3 fun betaq0 Qcutoff useindx
clear factor12 factor23 checkpoints Ncheckpoints Basis A1 b1 C1 d1 Aeq1 beq1 lb1 ub1 A2 b2 C2 d2 Aeq2 beq2 lb2 ub2 A3 b3 C3 d3 Aeq3 beq3 lb3 ub3 A b C d
 

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% NOW COMPUTE KNOT VECTOR AND SHAPE CONSTRAINTS FOR THE B-SPLINE COST FUNCTION:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isempty(knotvec_c)
    default = 7 ;

    K_c            = 7 ;  %%%%K_c=7 was chosen because further subintervals do not change point estimates in an economically important way.
    knotvec_c      = linspace(0,1,K_c+1)' ;
    [Gtemp,Ttemp]  = ecdf(AchievementData.T(AchievementData.Q>=1)) ;
    knotvec_c(1)   = 0 ;
    knotvec_c(end) = Ttemp(end) ;
    Gtemp(1)       = [] ;
    Ttemp(1)       = [] ;
    knotvec_c(2:end-1) = interp1(Gtemp,Ttemp,knotvec_c(2:end-1),'linear') ;
else
    default     = length(knotvec_c) - 1 ;  %#ok
    L.knotvec_c = knotvec_c ;
    K_c         = length(knotvec_c) - 1 ;
    L.K_c       = K_c ;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  Monotonicity Constraints:  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Aineq = [] ;
for ii=1:K_c+2
    Aineq = [Aineq;
        zeros(1,ii-1),1,-1,zeros(1,K_c+2-ii)]; %#ok
end
bineq = -tol*ones(K_c+2,1) ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  Convexity Constraints:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%THIS IS THE NEW WAY OF DOING IT, USING SIMPLE LINEAR CONSTRAINTS ON THE PARAMETERS:
Gamma = knotvec_c([(2:K_c+1)';K_c+1;K_c+1]) - knotvec_c([1;1;(1:K_c)']) ;
AineqConvexity = zeros(K_c+1,K_c+3) ;
bineqConvexity = zeros(K_c+1,1) ;
for jj=1:K_c+1
    AineqConvexity(jj,jj)   = -1/Gamma(jj) ;
    AineqConvexity(jj,jj+1) =  1/Gamma(jj) + 1/Gamma(jj+1) ;
    AineqConvexity(jj,jj+2) = -1/Gamma(jj+1) ;
    bineqConvexity(jj)      = -tol ;
end
Aineq = [Aineq; AineqConvexity] ;
bineq = [bineq; bineqConvexity] ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  Boundary Condition: 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%The following forces the value of the cost function at t=0 to be 0
Aeq = [1, zeros(1,K_c+2)] ;
beq = 0 ;
%%%%Normalization: The following forces the derivative of the cost function at t=0 to be 1

%%%%THE NEW WAY OF DOING IT, USING SIMPLE LINEAR CONSTRAINTS ON THE PARAMETERS
Aeq = [Aeq; 0, 1, zeros(1,K_c+1)] ; 
beq = [beq; (knotvec_c(2)-knotvec_c(1))/3] ; 


if ~isempty(initialguess)
    params0 = initialguess ;
else
    params0 = initialguessC(K_c) ;
end


lb         = -inf*ones(1,K_c+3) ; 
ub         = inf*ones(1,K_c+3) ; 
%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  FOR CONTINUATION-VALUE COMPUTATIONS:  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%For the domain of EQgivenT, we start with the knot vector, and then fill it in somewhat by adding midpoints.  
%%

Tdom     = sort([(2:2:16)'; knotvec_c]) ;

%%%%The array EQgivenT will contain expected future successes, given future time spent, t, when times are scaled by ThetaE.
%%%%The EQgivenT_ub version has 5 pages, and assumes ThetaE is the lower-bound (i.e., highest productivity) ThetaE-type
%%%%within each segment.  The EQgivenT_lb version assumes ThetaE is the upper-bound (i.e., lowest productivity) ThetaE-type
%%%%within each segment.
EQgivenT_ub        = nan(length(Tdom),Qbar+1,length(UdomHet(1,:))) ; %%%%THIS IS THE UPPER BOUND ON PAYOFF PER UNIT TIME (I.E., USING min(ThetaE)) WITHIN EACH SEGMENT.  The first column represents a baseline of zero completed tasks...
EQgivenT_ub(1,:,:) = 0 ; %%%With zero extra time, the continuation payoff is trivially zero.
EQgivenT_lb        = EQgivenT_ub ; %%%%THIS IS THE LOWER BOUND ON PAYOFF PER UNIT TIME (I.E., USING max(ThetaE)) WITHIN EACH SEGMENT.  The first column represents a baseline of zero completed tasks...
NSamps             = 1000 ; 

cutoffs = EPointEst.cutoffs ;
for seg = 1:length(UdomHet(1,:)) 
    F_utemp        = F_usmoothHet(:,seg) ; 
    Udomtemp       = UdomHet(:,seg) ; 
    flag           = find(abs(F_utemp(2:end)-F_utemp(1:end-1))<10^-6 | abs(Udomtemp(2:end)-Udomtemp(1:end-1))<10^-6) ; 
    F_utemp(flag)  = [] ; 
    Udomtemp(flag) = [] ; 
    tempUsamp = reshape(interp1(F_utemp,Udomtemp,rand(Qbar*NSamps,1),'linear'),Qbar,NSamps) ; 
    for aa=1:length(EQgivenT_ub(1,:,1)) 
        if aa==1
            tempUsamp2 = cutoffs(seg)*tau0*cumsum((tau1.^([1;zeros(Qbar-1,1)]*ones(1,NSamps))).*tempUsamp.*( ((aa:aa+Qbar-1).^-EPointEst.phi)'*ones(1,NSamps) )) ;
        else
            tempUsamp2 = cutoffs(seg)*tau0*cumsum(tempUsamp.*( ((aa:aa+Qbar-1).^-EPointEst.phi)'*ones(1,NSamps) )) ;
        end

        for tt=2:length(Tdom) 
            A = tempUsamp2 <= Tdom(tt) ; 
            [MaxVal,Amax] = max(flip(A)) ; 
            EQgivenT_ub(tt,aa,seg) = max(mean((Qbar+1-Amax).*(MaxVal>0)),Tdom(tt)*(NSamps*10)^-1) ;  
        end 
        tempUsamp2 = (cutoffs(seg+1)/cutoffs(seg))*tempUsamp2 ; 
        for tt=2:length(Tdom) 
            A = tempUsamp2 <= Tdom(tt) ; 
            [MaxVal,Amax] = max(flip(A)) ; 
            EQgivenT_lb(tt,aa,seg) = max(mean((Qbar+1-Amax).*(MaxVal>0)), Tdom(tt)*(NSamps*10)^-1) ;  
        end 
    end 
end 

if PreliminariesOnly==1
    if bootstrap(1)==0 || bootstrap(1)==10 || bootstrap(1)==11
        L.knotvec_c       = knotvec_c ;
        L.student_id      = student_id_l ;
        L.contract        = contract ;
        L.pay             = pay ;
        L.Q               = Q ;
        L.T               = T ;
        L.ThetaE          = ThetaE ;
        L.Activeids       = Activeids ;
        L.Marginalids1    = Marginalids1 ;
        L.Marginalids2    = Marginalids2 ;
        L.Q80freq         = Q80freq ;
        L.SimTimeCum      = SimTimeCum ;
        L.betaq1          = betaq1 ;
        L.betaq2          = betaq2 ;
        L.betaq3          = betaq3 ;
        L.knotvec_q       = knotvec_q ;
        L.Tdom            = Tdom ;
        L.EQgivenT_lb     = EQgivenT_lb ;
        L.EQgivenT_ub     = EQgivenT_ub ;
        L.cutoffs         = cutoffs ;
        L.Aeq             = Aeq ;
        L.beq             = beq ;
        L.Aineq           = Aineq ;
        L.bineq           = bineq ;
        L.lb              = lb ;
        L.ub              = ub ;
        
    elseif bootstrap(1)>=1
        L.knotvec_c       = knotvec_c ;
        L.student_id      = student_id_l ;
        L.contract        = contract ;
        L.pay             = pay ;
        L.Q               = Q ;
        L.T               = T ;
        L.ThetaE          = ThetaE ;
        L.Activeids       = Activeids ;
        L.Marginalids1    = Marginalids1 ;
        L.Marginalids2    = Marginalids2 ;
        if bootstrap(1)==4
            Q80freq(1,1,:)     = 80 ; 
            Q80freq(1,2,1)     = sum(Q80freq(:,2,1),'omitnan') ; 
            Q80freq(1,2,2)     = sum(Q80freq(:,2,2),'omitnan') ; 
            Q80freq(1,2,3)     = sum(Q80freq(:,2,3),'omitnan') ; 
            Q80freq(2:end,:,:) = nan ; 
        end
        L.Q80freq         = Q80freq ;
        L.SimTimeCum      = SimTimeCum ;
        L.Tdom            = Tdom ;
        L.EQgivenT_lb     = EQgivenT_lb ;
        L.EQgivenT_ub     = EQgivenT_ub ;
        L.cutoffs         = cutoffs ;
        L.Aeq             = Aeq ;
        L.beq             = beq ;
        L.Aineq           = Aineq ;
        L.bineq           = bineq ;
        L.lb              = lb ;
        L.ub              = ub ;
    end
else
    if bootstrap(1)==0 || bootstrap(1)==10 || bootstrap(1)==11
        myfun = @(x) ObjFun_cMSM(x,knotvec_c,student_id_l,Q,T,ThetaE,contract,Q80freq,pay,SimTimeCum,Activeids,Marginalids1,Marginalids2,...
            Tdom,EQgivenT_lb,EQgivenT_ub,cutoffs,betaq1,betaq2,betaq3,knotvec_q,0,bootstrap(1)) ;
    elseif bootstrap(1)>=1 && bootstrap(1)<10
        myfun = @(x) ObjFun_cMSM(x,knotvec_c,student_id_l,Q,T,ThetaE,contract,Q80freq,pay,SimTimeCum,Activeids,Marginalids1,Marginalids2,...
            Tdom,EQgivenT_lb,EQgivenT_ub,cutoffs,[],[],[],[],0,1) ;
    end
    %%%%HERE IS A LIST OF THE INPUTS TO THE OBJECTIVE FUNCTION:
    %%%% X (first input): this is the parameter vector for the cost function
    %%%% KNOTVEC_C: this is the knot vector for the B-spline Cost function
    %%%% STUDENT_ID_L, Q, T, CONTRACT: student-level variables, with rows of each corresponding to the same student.
    %%%% PAY:
    %%%% Q80FREQ: this is a 5x2x3 table.  Each sheet contains the extrapolated aggregate work values in column 1, and frequencies in column 2.  Sheet i is for contract group i.
    %%%% SIMTIMECUM: This contains the cumulative sum (within columns) of SimTime,  These are scaled by
    %%%% ThetaE, so the qth row represents a simulated cumulative work time through q units of success.
    %%%% Each column is a separate simulation, and each sheet is for an individual student.
    %%%% ACTIVEIDS: student id numbers for all who completed at least 2 learning tasks
    %%%% MARGINALIDS1 and MARGINALIDS2: student id numbers for all in contract groups 1 and 2 who completed exactly 1 learning task
    %%%% CUTOFFS: contains the cutoffs in ThetaE space for the heteroskedastic shock distributions
    %%%% DISPLAY (final input): determines whether plots are generated.
    %%%%
    %%%%
    %%%% SimTimeCumthreshold, SimTimeCum0, and SimNorm were used for the OUTDATED SELECTION CORRECTION
    %%%% G1emp, G2emp, G3emp, and Qemp were used for the OUTDATED empirical CDFs.  We now use the continuous versions given by betaq1, betaq2, betaq3, and knotvecq
    
    
    
    tol_solver = 10^-8 ;
    if bootstrap(1)==0
        
        disp('***')
        disp('***')
        disp('***')
        disp('Now estimating B-Spline Cost Function...')
        options_cpattern = optimoptions('patternsearch','InitialMeshSize',64,'UseCompleteSearch',true,... % 'MaxIterations',150, 'MaxFunEvals',2000,...
            'ConstraintTolerance',tol_solver, 'Display','iter', 'UseParallel',true,...
            'UseCompletePoll',true,'PlotFcn',@psplotbestf,'PollMethod','GSSPositiveBasis2N') %#ok
        [paramspat2,ssrvalcpat2,exitflagcpat2,outputcpat2] = patternsearch(myfun,params0,Aineq,bineq,Aeq,beq,lb,ub,[],options_cpattern) %#ok
        

        L.params          = paramspat2 ;
        L.knotvec_c       = knotvec_c ;
        L.exitflagcpat2   = exitflagcpat2 ;
        L.outputcpat2     = outputcpat2 ;
        L.ssrvalc         = ssrvalcpat2 ;
        
        [~,ObjFunOut]     = ObjFun_cMSM(paramspat2,knotvec_c,student_id_l,Q,T,ThetaE,contract,Q80freq,pay,SimTimeCum,Activeids,Marginalids1,Marginalids2,...
                                         Tdom,EQgivenT_lb,EQgivenT_ub,cutoffs,betaq1,betaq2,betaq3,knotvec_q,1,0) ;
        L.student_id      = ObjFunOut.student_id ;
        L.contract        = ObjFunOut.contract ;
        L.pay             = pay ;
        L.Q               = ObjFunOut.Q ;
        L.T               = ObjFunOut.T ;
        L.ThetaE          = ObjFunOut.ThetaE ;
        L.ssr_1lowtail    = ObjFunOut.ssr_1lowtail ;
        L.ssr_2lowtail    = ObjFunOut.ssr_2lowtail ;
        L.ssr_3lowtail    = ObjFunOut.ssr_3lowtail ;
        L.ssr_1mid        = ObjFunOut.ssr_1mid ;
        L.ssr_2mid        = ObjFunOut.ssr_2mid ;
        L.ssr_3mid        = ObjFunOut.ssr_3mid ;
        L.ssr_1upptail    = ObjFunOut.ssr_1upptail ;
        L.ssr_2upptail    = ObjFunOut.ssr_2upptail ;
        L.ssr_3upptail    = ObjFunOut.ssr_3upptail ;
        L.guardrail1up    = ObjFunOut.guardrail1up ;
        L.guardrail2up    = ObjFunOut.guardrail2up ;
        L.guardrail3up    = ObjFunOut.guardrail3up ;
        L.guardrail1lo    = ObjFunOut.guardrail1lo ;
        L.guardrail2lo    = ObjFunOut.guardrail2lo ;
        L.guardrail3lo    = ObjFunOut.guardrail3lo ;
        L.ssr             = ObjFunOut.ssr ;
        L.Activeids       = Activeids ;
        L.ThetaLACT       = ObjFunOut.ThetaLACT ;
        L.Marginalids1    = Marginalids1 ;
        L.ThetaLMRG1      = ObjFunOut.ThetaLMRG1 ;
        L.Marginalids2    = Marginalids2 ;
        L.ThetaLMRG2      = ObjFunOut.ThetaLMRG2 ;
        L.ThetaL80        = ObjFunOut.ThetaL80 ;
        L.Q80freq         = ObjFunOut.Q80freq ;
        L.SimTimeCum      = SimTimeCum ;
        L.betaq1          = betaq1 ;
        L.betaq2          = betaq2 ;
        L.betaq3          = betaq3 ;
        L.knotvec_q       = knotvec_q ;
        L.Tdom            = Tdom ;
        L.EQgivenT_lb     = EQgivenT_lb ;
        L.EQgivenT_ub     = EQgivenT_ub ;
        
    elseif bootstrap(1)>=1
        
        tol_solver = 10^-6 ;
        
        if bootstrap(1)==1
            %%%%FIRST PASS: USE A BIG INITIAL MESH SIZE, GSS algorithm, and allow incomplete polling
            options_cpattern = optimoptions('patternsearch','InitialMeshSize',4096,'UseCompleteSearch',true,... % 'MaxIterations',150, 'MaxFunEvals',2000,...
                'ConstraintTolerance',tol_solver,'FunctionTolerance',1e-4, 'MeshTolerance',1e-5,'StepTolerance',1e-5, 'Display','iter', 'UseParallel',true,...
                'PlotFcn',@psplotbestf,'PollMethod','GSSPositiveBasis2N') %#ok
            [paramspat2,ssrvalcpat2,exitflagcpat2,outputcpat2] = patternsearch(myfun,params0,Aineq,bineq,Aeq,beq,lb,ub,[],options_cpattern) ;
        elseif bootstrap(1)==2
            %%%%SECOND PASS: USE A SMALLER INITIAL MESH SIZE, GSS algorithm, and require complete polling
            options_cpattern = optimoptions('patternsearch','InitialMeshSize',16,'UseCompleteSearch',true,... % 'MaxIterations',150, 'MaxFunEvals',2000,...
                'ConstraintTolerance',tol_solver,'FunctionTolerance',1e-4, 'MeshTolerance',1e-5,'StepTolerance',1e-5,  'Display','iter', 'UseParallel',true,...
                'UseCompletePoll',true,'PlotFcn',@psplotbestf,'PollMethod','GSSPositiveBasis2N') %#ok
            [paramspat2,ssrvalcpat2,exitflagcpat2,outputcpat2] = patternsearch(myfun,params0,Aineq,bineq,Aeq,beq,lb,ub,[],options_cpattern) ;
        elseif bootstrap(1)==3
            %%%%THIRD PASS: USE A SMALLER INITIAL MESH SIZE, GPS algorithm, and require complete polling
            options_cpattern = optimoptions('patternsearch','InitialMeshSize',4,'UseCompleteSearch',true,... % 'MaxIterations',150, 'MaxFunEvals',2000,...
                'ConstraintTolerance',tol_solver,'FunctionTolerance',1e-4, 'MeshTolerance',1e-5,'StepTolerance',1e-5, 'Display','iter', 'UseParallel',true,...
                'UseCompletePoll',true,'PlotFcn',@psplotbestf,'PollMethod','GPSPositiveBasis2N') %#ok
            [paramspat2,ssrvalcpat2,exitflagcpat2,outputcpat2] = patternsearch(myfun,params0,Aineq,bineq,Aeq,beq,lb,ub,[],options_cpattern) ;
        
        
        elseif bootstrap(1)==10 %%%%THIS IS FOR COMPUTING POINT ESTIMATES WITH GENETIC ALGORITHM AND MAXGENERATIONS=15 
            options_cgenetic = optimoptions('ga','ConstraintTolerance',1e-6,'FunctionTolerance',1e-6,...
                'PlotFcn', {@gaplotbestf,@gaplotrange,@gaplotdistance},'MaxGenerations',20,'MaxStallGenerations',10,...
                'Display','iter','InitialPopulationMatrix',InitialPopMat,'PopulationSize',length(InitialPopMat(:,1)),'UseParallel',true) %#ok
            [paramspat2,ssrvalcpat2,exitflagcpat2,outputcpat2,population,scores] = ga(myfun,length(params0),Aineq,bineq,Aeq,beq,lb,ub,[],options_cgenetic) ;
        elseif bootstrap(1)==11 %%%%THIS IS FOR COMPUTING POINT ESTIMATES WITH GENETIC ALGORITHM AND MAXGENERATIONS=200
            options_cgenetic = optimoptions('ga','ConstraintTolerance',1e-6,'FunctionTolerance',1e-6,...
                'PlotFcn', {@gaplotbestf,@gaplotrange,@gaplotdistance},'MaxGenerations',200,'MaxStallGenerations',50,...
                'Display','iter','InitialPopulationMatrix',InitialPopMat,'PopulationSize',length(InitialPopMat(:,1)),'UseParallel',true) %#ok
            [paramspat2,ssrvalcpat2,exitflagcpat2,outputcpat2] = ga(myfun,length(params0),Aineq,bineq,Aeq,beq,lb,ub,[],options_cgenetic) ;
        
        
        elseif bootstrap(1)==12 %%%%THIS IS FOR COMPUTING BOOTSTRAPS WITH GENETIC ALGORITHM AND MAXGENERATIONS=15
            myfun = @(x) ObjFun_cMSM(x,knotvec_c,student_id_l,Q,T,ThetaE,contract,Q80freq,pay,SimTimeCum,Activeids,Marginalids1,Marginalids2,...
                Tdom,EQgivenT_lb,EQgivenT_ub,cutoffs,[],[],[],[],0,bootstrap(1)) ; 
            options_cgenetic = optimoptions('ga','ConstraintTolerance',1e-6,'FunctionTolerance',1e-4,...
                'PlotFcn', {@gaplotbestf,@gaplotrange,@gaplotdistance},'MaxGenerations',10,'MaxStallGenerations',5,...
                'Display','iter','InitialPopulationMatrix',InitialPopMat,'PopulationSize',length(InitialPopMat(:,1)),'UseParallel',true) %#ok
            [paramspat2,ssrvalcpat2,exitflagcpat2,outputcpat2,population,scores] = ga(myfun,length(params0),Aineq,bineq,Aeq,beq,lb,ub,[],options_cgenetic) ;
        elseif bootstrap(1)==13  %%%%THIS IS FOR COMPUTING BOOTSTRAPS WITH GENETIC ALGORITHM AND MAXGENERATIONS=100
            myfun = @(x) ObjFun_cMSM(x,knotvec_c,student_id_l,Q,T,ThetaE,contract,Q80freq,pay,SimTimeCum,Activeids,Marginalids1,Marginalids2,...
                                      Tdom,EQgivenT_lb,EQgivenT_ub,cutoffs,[],[],[],[],0,bootstrap(1)) ; 
            options_cgenetic = optimoptions('ga','ConstraintTolerance',1e-6,'FunctionTolerance',1e-4,...
                'PlotFcn', {@gaplotbestf,@gaplotrange,@gaplotdistance},'MaxGenerations',100,'MaxStallGenerations',15,...
                'Display','iter','InitialPopulationMatrix',InitialPopMat,'PopulationSize',length(InitialPopMat(:,1)),'UseParallel',true) %#ok 
            [paramspat2,ssrvalcpat2,exitflagcpat2,outputcpat2] = ga(myfun,length(params0),Aineq,bineq,Aeq,beq,lb,ub,[],options_cgenetic) ; 
        end
        
        
        L.params       = paramspat2 ;
        L.knotvec_c    = knotvec_c ;
        L.exitflagc    = exitflagcpat2 ;
        L.outputcpat2  = outputcpat2 ;
        L.ssrvalc      = ssrvalcpat2 ;
        if bootstrap(1)==10 || bootstrap(1)==12
            L.population = population ;
            L.scores     = scores ;
        end
        
        if bootstrap(1)==0 || bootstrap(1)==10 || bootstrap(1)==11 
            [~,ObjFunOut]  = ObjFun_cMSM(paramspat2,knotvec_c,student_id_l,Q,T,ThetaE,contract,...
                                          Q80freq,pay,SimTimeCum,Activeids,Marginalids1,Marginalids2,...
                                          Tdom,EQgivenT_lb,EQgivenT_ub,cutoffs,betaq1,betaq2,betaq3,knotvec_q,0,bootstrap(1)) ; 
        else
            [~,ObjFunOut]  = ObjFun_cMSM(paramspat2,knotvec_c,student_id_l,Q,T,ThetaE,contract,...
                                          Q80freq,pay,SimTimeCum,Activeids,Marginalids1,Marginalids2,...
                                          Tdom,EQgivenT_lb,EQgivenT_ub,cutoffs,[],[],[],[],0,bootstrap(1)) ; 
        end
        L.student_id      = ObjFunOut.student_id ; 
        L.contract        = ObjFunOut.contract ; 
        L.pay             = pay ; 
        L.Q               = ObjFunOut.Q ; 
        L.T               = ObjFunOut.T ; 
        L.ThetaE          = ObjFunOut.ThetaE ; 
        L.ssr_1lowtail    = ObjFunOut.ssr_1lowtail ; 
        L.ssr_2lowtail    = ObjFunOut.ssr_2lowtail ; 
        L.ssr_3lowtail    = ObjFunOut.ssr_3lowtail ; 
        L.ssr_1mid        = ObjFunOut.ssr_1mid ; 
        L.ssr_2mid        = ObjFunOut.ssr_2mid ; 
        L.ssr_3mid        = ObjFunOut.ssr_3mid ; 
        L.ssr_1upptail    = ObjFunOut.ssr_1upptail ; 
        L.ssr_2upptail    = ObjFunOut.ssr_2upptail ; 
        L.ssr_3upptail    = ObjFunOut.ssr_3upptail ; 
        L.guardrail1up    = ObjFunOut.guardrail1up ; 
        L.guardrail2up    = ObjFunOut.guardrail2up ; 
        L.guardrail3up    = ObjFunOut.guardrail3up ; 
        L.guardrail1lo    = ObjFunOut.guardrail1lo ; 
        L.guardrail2lo    = ObjFunOut.guardrail2lo ; 
        L.guardrail3lo    = ObjFunOut.guardrail3lo ; 
        L.ssr             = ObjFunOut.ssr ; 
        L.Activeids       = Activeids ; 
        L.ThetaLACT       = ObjFunOut.ThetaLACT ; 
        L.Marginalids1    = Marginalids1 ; 
        L.ThetaLMRG1      = ObjFunOut.ThetaLMRG1 ; 
        L.Marginalids2    = Marginalids2 ; 
        L.ThetaLMRG2      = ObjFunOut.ThetaLMRG2 ; 
        L.ThetaL80        = ObjFunOut.ThetaL80 ; 
        L.Q80freq         = ObjFunOut.Q80freq ; 
        L.SimTimeCum      = SimTimeCum ; 
        L.Tdom            = Tdom ; 
        L.EQgivenT_lb     = EQgivenT_lb ; 
        L.EQgivenT_ub     = EQgivenT_ub ; 
    end
    
    
    
    

    
    
    
    
    
end
end





